Excercise 4

access the MASS package

library(MASS)

load the data

data("Boston")

explore the dataset

str(Boston)

## 'data.frame':    506 obs. of  14 variables:
##  $ crim   : num  0.00632 0.02731 0.02729 0.03237 0.06905 ...
##  $ zn     : num  18 0 0 0 0 0 12.5 12.5 12.5 12.5 ...
##  $ indus  : num  2.31 7.07 7.07 2.18 2.18 2.18 7.87 7.87 7.87 7.87 ...
##  $ chas   : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ nox    : num  0.538 0.469 0.469 0.458 0.458 0.458 0.524 0.524 0.524 0.524 ...
##  $ rm     : num  6.58 6.42 7.18 7 7.15 ...
##  $ age    : num  65.2 78.9 61.1 45.8 54.2 58.7 66.6 96.1 100 85.9 ...
##  $ dis    : num  4.09 4.97 4.97 6.06 6.06 ...
##  $ rad    : int  1 2 2 3 3 3 5 5 5 5 ...
##  $ tax    : num  296 242 242 222 222 222 311 311 311 311 ...
##  $ ptratio: num  15.3 17.8 17.8 18.7 18.7 18.7 15.2 15.2 15.2 15.2 ...
##  $ black  : num  397 397 393 395 397 ...
##  $ lstat  : num  4.98 9.14 4.03 2.94 5.33 ...
##  $ medv   : num  24 21.6 34.7 33.4 36.2 28.7 22.9 27.1 16.5 18.9 ...

summary(Boston)

##       crim                zn             indus            chas        
##  Min.   : 0.00632   Min.   :  0.00   Min.   : 0.46   Min.   :0.00000  
##  1st Qu.: 0.08204   1st Qu.:  0.00   1st Qu.: 5.19   1st Qu.:0.00000  
##  Median : 0.25651   Median :  0.00   Median : 9.69   Median :0.00000  
##  Mean   : 3.61352   Mean   : 11.36   Mean   :11.14   Mean   :0.06917  
##  3rd Qu.: 3.67708   3rd Qu.: 12.50   3rd Qu.:18.10   3rd Qu.:0.00000  
##  Max.   :88.97620   Max.   :100.00   Max.   :27.74   Max.   :1.00000  
##       nox               rm             age              dis        
##  Min.   :0.3850   Min.   :3.561   Min.   :  2.90   Min.   : 1.130  
##  1st Qu.:0.4490   1st Qu.:5.886   1st Qu.: 45.02   1st Qu.: 2.100  
##  Median :0.5380   Median :6.208   Median : 77.50   Median : 3.207  
##  Mean   :0.5547   Mean   :6.285   Mean   : 68.57   Mean   : 3.795  
##  3rd Qu.:0.6240   3rd Qu.:6.623   3rd Qu.: 94.08   3rd Qu.: 5.188  
##  Max.   :0.8710   Max.   :8.780   Max.   :100.00   Max.   :12.127  
##       rad              tax           ptratio          black       
##  Min.   : 1.000   Min.   :187.0   Min.   :12.60   Min.   :  0.32  
##  1st Qu.: 4.000   1st Qu.:279.0   1st Qu.:17.40   1st Qu.:375.38  
##  Median : 5.000   Median :330.0   Median :19.05   Median :391.44  
##  Mean   : 9.549   Mean   :408.2   Mean   :18.46   Mean   :356.67  
##  3rd Qu.:24.000   3rd Qu.:666.0   3rd Qu.:20.20   3rd Qu.:396.23  
##  Max.   :24.000   Max.   :711.0   Max.   :22.00   Max.   :396.90  
##      lstat            medv      
##  Min.   : 1.73   Min.   : 5.00  
##  1st Qu.: 6.95   1st Qu.:17.02  
##  Median :11.36   Median :21.20  
##  Mean   :12.65   Mean   :22.53  
##  3rd Qu.:16.95   3rd Qu.:25.00  
##  Max.   :37.97   Max.   :50.00

library(tidyr)
library(corrplot)

## corrplot 0.84 loaded

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following object is masked from 'package:MASS':
## 
##     select

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

calculate the correlation matrix and round it

cor_matrix<-cor(Boston) %>% round(digits = 2 )

print the correlation matrix

cor_matrix

##          crim    zn indus  chas   nox    rm   age   dis   rad   tax
## crim     1.00 -0.20  0.41 -0.06  0.42 -0.22  0.35 -0.38  0.63  0.58
## zn      -0.20  1.00 -0.53 -0.04 -0.52  0.31 -0.57  0.66 -0.31 -0.31
## indus    0.41 -0.53  1.00  0.06  0.76 -0.39  0.64 -0.71  0.60  0.72
## chas    -0.06 -0.04  0.06  1.00  0.09  0.09  0.09 -0.10 -0.01 -0.04
## nox      0.42 -0.52  0.76  0.09  1.00 -0.30  0.73 -0.77  0.61  0.67
## rm      -0.22  0.31 -0.39  0.09 -0.30  1.00 -0.24  0.21 -0.21 -0.29
## age      0.35 -0.57  0.64  0.09  0.73 -0.24  1.00 -0.75  0.46  0.51
## dis     -0.38  0.66 -0.71 -0.10 -0.77  0.21 -0.75  1.00 -0.49 -0.53
## rad      0.63 -0.31  0.60 -0.01  0.61 -0.21  0.46 -0.49  1.00  0.91
## tax      0.58 -0.31  0.72 -0.04  0.67 -0.29  0.51 -0.53  0.91  1.00
## ptratio  0.29 -0.39  0.38 -0.12  0.19 -0.36  0.26 -0.23  0.46  0.46
## black   -0.39  0.18 -0.36  0.05 -0.38  0.13 -0.27  0.29 -0.44 -0.44
## lstat    0.46 -0.41  0.60 -0.05  0.59 -0.61  0.60 -0.50  0.49  0.54
## medv    -0.39  0.36 -0.48  0.18 -0.43  0.70 -0.38  0.25 -0.38 -0.47
##         ptratio black lstat  medv
## crim       0.29 -0.39  0.46 -0.39
## zn        -0.39  0.18 -0.41  0.36
## indus      0.38 -0.36  0.60 -0.48
## chas      -0.12  0.05 -0.05  0.18
## nox        0.19 -0.38  0.59 -0.43
## rm        -0.36  0.13 -0.61  0.70
## age        0.26 -0.27  0.60 -0.38
## dis       -0.23  0.29 -0.50  0.25
## rad        0.46 -0.44  0.49 -0.38
## tax        0.46 -0.44  0.54 -0.47
## ptratio    1.00 -0.18  0.37 -0.51
## black     -0.18  1.00 -0.37  0.33
## lstat      0.37 -0.37  1.00 -0.74
## medv      -0.51  0.33 -0.74  1.00

visualize the correlation matrix

corrplot(cor_matrix, method="circle", type="upper", cl.pos="b", tl.pos="d", tl.cex = 0.6)

It seems that “dis” is closely related to “age”, “nox”, “indus” and “zn”.

center and standardize variables

boston_scaled <- scale(Boston)

summaries of the scaled variables

summary(boston_scaled)

##       crim                 zn               indus        
##  Min.   :-0.419367   Min.   :-0.48724   Min.   :-1.5563  
##  1st Qu.:-0.410563   1st Qu.:-0.48724   1st Qu.:-0.8668  
##  Median :-0.390280   Median :-0.48724   Median :-0.2109  
##  Mean   : 0.000000   Mean   : 0.00000   Mean   : 0.0000  
##  3rd Qu.: 0.007389   3rd Qu.: 0.04872   3rd Qu.: 1.0150  
##  Max.   : 9.924110   Max.   : 3.80047   Max.   : 2.4202  
##       chas              nox                rm               age         
##  Min.   :-0.2723   Min.   :-1.4644   Min.   :-3.8764   Min.   :-2.3331  
##  1st Qu.:-0.2723   1st Qu.:-0.9121   1st Qu.:-0.5681   1st Qu.:-0.8366  
##  Median :-0.2723   Median :-0.1441   Median :-0.1084   Median : 0.3171  
##  Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.0000  
##  3rd Qu.:-0.2723   3rd Qu.: 0.5981   3rd Qu.: 0.4823   3rd Qu.: 0.9059  
##  Max.   : 3.6648   Max.   : 2.7296   Max.   : 3.5515   Max.   : 1.1164  
##       dis               rad               tax             ptratio       
##  Min.   :-1.2658   Min.   :-0.9819   Min.   :-1.3127   Min.   :-2.7047  
##  1st Qu.:-0.8049   1st Qu.:-0.6373   1st Qu.:-0.7668   1st Qu.:-0.4876  
##  Median :-0.2790   Median :-0.5225   Median :-0.4642   Median : 0.2746  
##  Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.0000  
##  3rd Qu.: 0.6617   3rd Qu.: 1.6596   3rd Qu.: 1.5294   3rd Qu.: 0.8058  
##  Max.   : 3.9566   Max.   : 1.6596   Max.   : 1.7964   Max.   : 1.6372  
##      black             lstat              medv        
##  Min.   :-3.9033   Min.   :-1.5296   Min.   :-1.9063  
##  1st Qu.: 0.2049   1st Qu.:-0.7986   1st Qu.:-0.5989  
##  Median : 0.3808   Median :-0.1811   Median :-0.1449  
##  Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.0000  
##  3rd Qu.: 0.4332   3rd Qu.: 0.6024   3rd Qu.: 0.2683  
##  Max.   : 0.4406   Max.   : 3.5453   Max.   : 2.9865

class of the boston_scaled object

class(boston_scaled)

## [1] "matrix"

change the object to data frame

boston_scaled <- as.data.frame(boston_scaled)

summary of the scaled crime rate

summary(boston_scaled$crim)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## -0.419367 -0.410563 -0.390280  0.000000  0.007389  9.924110

create a quantile vector of crim and print it

bins <- quantile(boston_scaled$crim)

bins

##           0%          25%          50%          75%         100% 
## -0.419366929 -0.410563278 -0.390280295  0.007389247  9.924109610

create a categorical variable ‘crime’

crime <- cut(boston_scaled$crim, breaks = bins, include.lowest = TRUE, labels = c("low", "med_low", "med_high", "high"))

look at the table of the new factor crime

table(crime)

## crime
##      low  med_low med_high     high 
##      127      126      126      127

remove original crim from the dataset

boston_scaled <- dplyr::select(boston_scaled, -crim)

add the new categorical value to scaled data

boston_scaled <- data.frame(boston_scaled, crime)

number of rows in the Boston dataset

n <- nrow(Boston)

choose randomly 80% of the rows

ind <- sample(n,  size = n * 0.8)

create train set

train <- boston_scaled[ind,]

create test set

test <- boston_scaled[-ind,]

save the correct classes from test data

correct_classes <- test$crime

remove the crime variable from test data

test <- dplyr::select(test, -crime)

linear discriminant analysis

lda.fit <- lda(crime ~ ., data = train)

print the lda.fit object

lda.fit

## Call:
## lda(crime ~ ., data = train)
## 
## Prior probabilities of groups:
##       low   med_low  med_high      high 
## 0.2450495 0.2524752 0.2549505 0.2475248 
## 
## Group means:
##                  zn      indus        chas        nox          rm
## low       1.0461588 -0.9193671 -0.15302300 -0.8658565  0.47114176
## med_low  -0.1122754 -0.2836581 -0.04073494 -0.5616129 -0.12080475
## med_high -0.3939928  0.1320544  0.18636222  0.3577097  0.08906587
## high     -0.4872402  1.0171519 -0.11484506  1.0507196 -0.48334202
##                 age        dis        rad        tax     ptratio
## low      -0.9127575  0.8930324 -0.6860539 -0.7486573 -0.48988951
## med_low  -0.3200213  0.3202680 -0.5562628 -0.4490306 -0.03787817
## med_high  0.3436602 -0.3359473 -0.3797623 -0.3155899 -0.26557300
## high      0.8177631 -0.8520027  1.6377820  1.5138081  0.78037363
##               black       lstat        medv
## low       0.3843834 -0.77397494  0.55776469
## med_low   0.3197696 -0.14842695 -0.01305424
## med_high  0.0792625  0.01082824  0.17800131
## high     -0.7442657  0.90407242 -0.73314451
## 
## Coefficients of linear discriminants:
##                  LD1         LD2         LD3
## zn       0.077141148  0.89285925 -0.76600497
## indus   -0.003590024 -0.20471096  0.17344390
## chas    -0.080282919 -0.14037266  0.04955767
## nox      0.371905492 -0.52724127 -1.51518917
## rm      -0.076869255 -0.07326862 -0.15838327
## age      0.266279743 -0.23736386 -0.04097211
## dis     -0.059161664 -0.27421668  0.01015961
## rad      3.028258006  0.88930967 -0.28923422
## tax     -0.059480822 -0.07890009  0.93744527
## ptratio  0.112331132  0.03454140 -0.24823622
## black   -0.123100001  0.01009864  0.12053335
## lstat    0.221916198 -0.28380724  0.12044056
## medv     0.144970969 -0.35323617 -0.36984297
## 
## Proportion of trace:
##    LD1    LD2    LD3 
## 0.9428 0.0415 0.0156

the function for lda biplot arrows

lda.arrows <- function(x, myscale = 1, arrow_heads = 0.1, color = "orange", tex = 0.75, choices = c(1,2)){
  heads <- coef(x)
  arrows(x0 = 0, y0 = 0, 
         x1 = myscale * heads[,choices[1]], 
         y1 = myscale * heads[,choices[2]], col=color, length = arrow_heads)
  text(myscale * heads[,choices], labels = row.names(heads), 
       cex = tex, col=color, pos=3)
}

target classes as numeric

classes <- as.numeric(train$crime)

plot the lda results

plot(lda.fit, dimen = 2, col = classes, pch = classes)
lda.arrows(lda.fit, myscale = 1)

predict classes with test data

lda.pred <- predict(lda.fit, newdata = test)

cross tabulate the results

table(correct = correct_classes, predicted = lda.pred$class)

##           predicted
## correct    low med_low med_high high
##   low       13      15        0    0
##   med_low    5      15        4    0
##   med_high   0       3       20    0
##   high       0       0        0   27

load MASS and Boston

library(MASS)
data('Boston')

center and standardize variables

boston_scaled <- scale(Boston)

change the object to data frame from matrix type.

boston_scaled <- as.data.frame(boston_scaled)

Calculate the Euclidean distances between observations.

dist_eu <- dist(boston_scaled)

look at the summary of the distances

summary(dist_eu)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.1343  3.4625  4.8241  4.9111  6.1863 14.3970

K-means clustering

km <-kmeans(boston_scaled, centers = 3)

plot the Boston dataset with clusters

pairs(boston_scaled, col = km$cluster)

Boston dataset is available

set.seed(123)

determine the number of clusters

k_max <- 10

calculate the total within sum of squares

twcss <- sapply(1:k_max, function(k){kmeans(boston_scaled, k)$tot.withinss})

visualize the results

library(ggplot2)

qplot(x = 1:k_max, y = twcss, geom = 'line')

k-means clustering

km <-kmeans(boston_scaled, centers = 2)

plot the Boston dataset with clusters

pairs(boston_scaled, col = km$cluster)

ggpairs

library(GGally)

## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2

## 
## Attaching package: 'GGally'

## The following object is masked from 'package:dplyr':
## 
##     nasa

ggpairs(boston_scaled)

Super bonus

model_predictors <- dplyr::select(train, -crime)

check the dimensions

dim(model_predictors)

## [1] 404  13

dim(lda.fit$scaling)

## [1] 13  3

matrix multiplication

matrix_product <- as.matrix(model_predictors) %*% lda.fit$scaling
matrix_product <- as.data.frame(matrix_product)

library(plotly)

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:MASS':
## 
##     select

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

3D figure

plot_ly(x = matrix_product$LD1, y = matrix_product$LD2, z = matrix_product$LD3, type= 'scatter3d', mode='markers', color = train$crime )

First one needs to do k-means with 4 clusters to compare the methods.

km3D <-kmeans(boston_scaled, centers = 4)

3D figure

plot_ly(x = matrix_product$LD1, y = matrix_product$LD2, z = matrix_product$LD3, type= 'scatter3d', mode='markers', color = km3D$cluster[ind])